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Abstract 

Based on a geometric discretization scheme for Maxwell equations, we unveil a mathematical 
transformation between the electric field intensity E and the magnetic field intensity H, denoted 
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second-order vector wave equations. We show that the primal formulation recovers the conventional 
(edge-element) finite element method (FEM) and suggests a geometric foundation for it. On the 
other hand, the dual formulation suggests a new (dual) type of FEM. Although both formulations 
give identical dynamical physical solutions, the dimensions of the null spaces are different. 
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I. INTRODUCTION 



The finite element method (FEM), originally developed for structure design and analysis, 
is usually based on nodal elements Simply applying nodal elements to Maxwell equations 
causes problems such as spurious modes [2| . The use of edge elements is the only reasonable 
way |3[ to remove the spurious modes because the electric field intensity E is a differential 
1-form with degrees of freedom (DoFs) associated with the edges of a lattice 1 . 

The basic strategy of traditional FEM (Galerkin's method) is to seek the solution by 
weighting the residual of the second-order wave equations. Here, we adopt a different route. 
Based on a general discretization scheme for Maxwell equations on irregular lattices, we 
construct two system matrices in terms of the electric field intensity E (denoted as primal 
formulation) and the magnetic field intensity H (denoted as dual formulation), respectively. 
The primal formulation recovers the FEM based on edge elements, and suggests a geometric 
foundation for it. On the other hand, the dual formulation suggests a new (dual) type 
of FEM. Although both formulations give identical physical solutions, the dimensions of 
the null spaces are different. The connection between the primal formulation and dual 
formulation is established via a transformation denoted here as Galerkin duality (not to be 
confused with conventional electromagnetic duality 00). 



II. DISCRETE MAXWELL EQUATIONS 

Maxwell equations in source-free, three-dimensional (3D) space (in the Fourier domain) 
are written in terms of differential forms [t| 3| as 

dE = iuB, dB = 0, dH = -iuD, dD = 0, (1) 

where E and H are electric and magnetic field intensity 1-forms, D and B are electric 
and magnetic flux 2-forms, and d is the (metric-free) exterior derivative operator. We use 
the convention e~~ lujt throughout this paper. Constitutive equations, which include all metric 
information, are written in terms of Hodge (star) operators (that fix an isomorphism between 
p-forms and (3 — p)-forms) 



1 For high order 1-forms |4j, DoFs of 1-forms could also associate with the faces and volumes, but do not 
associate with the nodes. Pointed out by one reviewer, recent work by Rapetti and Bossavit suggests 
that DoFs for high order 1-forms are still fundamentally associated with segments (small edges). 
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D = -k e E , H = -kp-iB. (2) 

By applying basic tools of algebraic topology and a discrete analog of differential forms, 
discrete electromagnetic theory can be constructed from first principles on a general (irreg- 
ular) primal/dual lattice (oriented cell-complex) pf. The discrete Maxwell equations read 
as 9| 

[d mrl }E=icoM, [d div }M=0, [d* curl ] W=-iuB, KJD=0, (3) 

where E, B, H, D are arrays of DoFs and [d cur i] , [ddi V ] , [d* ur t ] , [d* div ] are incidence ma- 
trices that encode the discrete exterior derivatives (discrete counterparts to the curl and 
divergence operators, distilled from their metric structure) on the primal and dual lattice, 
respectively. Due to the absence of metric structure, entries of the incidence matrices assume 
only {—1, 0, 1} values [3]. 

The discrete Hodge operators can be, in general, written as follows 

D = [* 6 ] E, H = B. (4) 

One approach to construct the Hodge matrices [* e ] and [*^-i] will be discussed in next Sec- 
tion. The Hodge matrices should be positive definite because Hodge operators are positive 
definite (in a Riemannian manifold). 



III. DISCRETE HODGE OPERATORS 

Let Q be a n-dimensional differentiable manifold and F p (Q) the space of forms of p- 
degree defined on it. If Q is endowed with a metric, then the Hodge operator * : rj — > ^ = ~ki] 
HQ is defined as a map of n e F p (O) to ^ G F n ~ p (Q) such that for any ijj e F p (D) 



$ A i = / i/; A *tj. (5) 

The Hodge operator defines (through a metric) an infinite dimensional inner product, de- 
noted as (if), rj) 

= / i> ^*v- (6) 
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For some form if) we can also define the Hodge square of ip by 

(if), if)) = [ i)A*if), (7) 



which is positive when the metric is positive definite. By applying (JJJ) to electric field and 
magnetic field, one can obtain constitutive relations in terms of Hodge operators in 3D 
Euclidean space R 3 as 

(E, E) = [ EAD= [ E A * e E, (8) 

JR? JR? 



Whitney forms 



[B,B) = / BAH= SAV-B. (9) 
jr? Jr? 



12| are the basic interpolants for discrete differential forms of various 
degrees defined over tetrahedra. Whitney forms can be expressed in term of the barycentric 
coordinates associated with each tetrahedron nodes (Q, Q, (k,(r) as 



™° = Ct, (io) 

< = &dQ ~ Qdd, (11) 

= 2 (ddQ A dC k + Qd( k A dd + (kdd A dQ , (12) 

% n I CidQ A d( k A rfCr - ( r d(i A dQ A d( k 

W ijkr — 6 

y +CfcG?Cr A A dCj - Cj^Cfc A rf Cr A dd 

(See the appendix for the basis functions over cubes). Accordingly, we use Whitney 1-forms 
as the interpolants for electric field intensity 1-form E, and Whitney 2-forms as interpolants 
for the magnetic flux 2-form B 

E = J2 e i,3 W lp B = J2 b ^ w li,k- ( 14 ) 

Note that the above expansions guarantee tangential continuity of E and normal continuity 
of B simultaneously. 

Using these basis functions and the Euclidean metric, matrix representations for the 
Hodge operators * e and can be constructed by combining Eq. (JSJ), Eq. Q and Eq. 
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In the above, matrix entry [*e]{(j.,) (i j)} comes from edge and edge (i>3ji an d matrix 
entry [Vr^-i]^^ (J'jV)} comes from face (i,j,k) and face \i,j,kj. These matrices denoted 



as Galerkin's discrete Hodges 



16(, or simply Galerkin's Hodges. 



IV. PRIMAL AND DUAL DISCRETE WAVE EQUATIONS 
A. Discrete wave equations 

From Eqs.(JBJ), (jlj) and ffTHfl . two discrete, second-order vector wave equations can be 
obtained 

[d* curl ][^-i][d curl ]E = u 2 [* e ]E, (16) 
[d C un] [* e ] _1 [d* curl ] M = u? [v,]" 1 H, (17) 

corresponding to a primal and dual formulation, respectively. These are the discrete analogs 
of the curl curl equations 

V*- x ^ x 2? = u 2 e~fi, (18) 
A* 

^- x ^ x ^ = u 2 fJf. (19) 
e 

It can be shown that [d* url ] [* M -i] [d CU ri} is identical to the conventional stiffness matrix [S] 
(see Appendix), arising in FEM using edge elements 

W{M,m = j I x ■ >< ^j) dv - ( 2 °) 

Moreover, the Hodge matrix [* e ] is identical to the conventional mass matrix. Hence, the 
primal formulation recovers the conventional edge-element FEM and suggests a geometric 
foundation for it. For the dual formulation, we can introduce dual stiffness [S^] and mass 
[M^l matrices 

[st] = MKr'tCi], (2i) 

[Mt] = fv-x]- 1 . (22) 

This dual formulation has no direct counterpart in traditional FEM. As discussed next, these 
two formulations lead to the same dynamic solutions, but have very different mathematical 
properties. Note that the Hodge matrices [* e ] and are sparse matrices, but their 

inverses [* e ] -1 and are in general not sparse. 



TABLE I: Galerkin duality vs. Electromagnetic duality. 



Galerkin duality 


Electromagnetic duality 






PEC^PEC 


PEC^PMC 


Dirichlet BC— ^Neumann BC 


Dirichlet BC^Dirichlet BC 


Neumann BC— ^Dirichlet BC 


Neumann BC^Neumann BC 



B. Galerkin duality 



Galerkin duality is a mathematical transformation between the above primal and dual 
formulations. Note that Galerkin duality is distinct from usual electromagnetic duality 
as illustrated in Table I. 



Based on Galerkin duality and the discrete Hodge operators introduced before, we can 
construct two different system matrices for a given problem 

[X E ] = i^]- 1 [d* mrl ] [d curl ] , (23) 

[X H ] = K-A [dcurl] [*J _1 • ( 24 ) 

Both [Xe] and [Xh] encode all discrete dynamic information, and hence produce identical 
dynamic solutions. However, their null spaces (associated with zero modes) are very differ- 
ent. In other words, for a discretization of the same physical system, the dimensions of the 
(discrete) zero eigenspaces are different under Galerkin duality. This can be explained by 
algebraic properties of discrete Hodge decomposition, and verified by numerical simulations, 
as discussed in Section V. 



C. An approach to handle Neumann boundary condition 

Since Dirichlet boundary condition and Neumann boundary condition are Galerkin dual 
to each other for some underlying differential equations, we propose an approach to handle 
Neumann boundary condition. Consider a differential equation 

90 = 0, (25) 



6 



where G is a differential operator and <p is the unknown physical quantity, with Neumann 
boundary condition. By Galerkin duality, this problem is equivalent to solving 



with Dirichlet boundary condition. Here 0t is the Galerkin dual to B, and is the Galerkin 
dual to 4>. Note that by using Galerkin duality, we can transform Eq. into Eq. (pjj) . since 
it may be much easier to handle Dirichlet boundary conditions than Neumann boundary 
conditions in some problems. The function <ft can be expanded in terms of basis functions 
W- (e.g., Whitney forms) as 



V. EXAMPLES 

To demonstrate the Galerkin duality, we provide some numerical simulations for 2D 
cavity problems in the (x, y) plane. Both TE and TM cases are simulated. The finite 
element meshes for these examples were generated by using Triangle, a freely available 2D 
mesh generator . The angular frequencies of the resonant modes are obtained by solving 
the eigenvalue equation (fTB|) (primal formulation) or the eigenvalue equation (JT7J) (dual 
formulation). For simplicity, we set e = \i = 1. 

A. Whitney forms in 2D 

The vector proxies of Whitney forms in 2D can be written in term of barycentric coordi- 
nates (&,(jXk) as 



eV f = o 



(26) 




(27) 




2 (ovo x vc* + ova x vc< + CfcVO x vo) . 



(28) 
(29) 
(30) 



In the above, W® and Wf, k are scalars and is a vector. 
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TABLE II: TE vs.TM. 





E 


B 


Degree of differential-form (TE) 


1 


2 


Degree of differential-form (TM) 





1 


Element (TE) 


edge 


face 


Element (TM) 


node 


edge 



1. TE case 

For the TE case, we use ■ as the interpolants for the electric field intensity ~E* and 
W?j k as the interpolants for the magnetic flux B z 

^ = E e M'C I> £ b id ,kW^ k . (31) 

Galerkin's Hodges become 

2. TM case 

For the TM case, we use Wf as the interpolants for the electric field intensity E z and 
W id as the interpolants for the magnetic flux ~B* 

Galerkin's Hodges become 

[*J{i,T} = / eW?-W?dS, 
^-] {(M ),(IJ)} = f^lr^hdS. (34) 
The comparison between TE and TM case is illustrated in Table II 
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TABLE III: TE modes (the angular frequencies of the 5 lowest nonzero modes) of a circular cavity. 



Mode TEmr, 


Primal 


Dual 


Analytical 


Error (%) 




1.8493 


1.8493 


1.8412 


0.4416 




1.8494 


1.8494 


1.8412 


0.4483 


TE%\ 


3.0707 


3.0707 


3.0542 


0.5381 


TE21 


3.0708 


3.0708 


3.0542 


0.5412 


TEqi 


3.8421 


3.8421 


3.8317 


0.2705 


# zero modes 


136 


1 






# nonzero modes 


311 


311 







B. Circular cavity 

Table III and Table IV present the results for TE modes and TM modes of a circular 
cavity with radius a — 1. The analytical solutions of TE modes are the zeros of Bessel 
function derivative J' m (x) ; The analytical solutions of TM modes are the zeros of Bessel 
function J m (x). Note that TE mn and TM mn have a twofold degeneracy analytically if 
m 7^ 0. However, the numerical solutions break the degeneracy. From the Table III (2D 
TE modes), we find that the number of zero modes of primal formulation is equal to the 
number of internal nodes, while the number of zero modes of dual formulation is 1. On the 
other hand, from the Table IV (2D TM modes), we find that the number of zero modes of 
primal formulation is 0, while the number of zero modes of dual formulation is Np — 1. From 
the last rows of Table III and Table IV, we conclude that both formulations give the same 
number of nonzero modes. These numerical facts, summarized in Table V, will be explained 
by applying a discrete Hodge decomposition in next subsection. 

C. Discrete Hodge decomposition 



as 



n a contractible domain ft, the Hodge decomposition for a p-form F p (ft) can be written 

F p (ft) = dFP- 1 (ft) © 5F P+1 (ft) , (35) 
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TABLE IV: TM modes (the angular frequencies of the 5 lowest nonzero modes) of a circular cavity. 



Mode TM mn 


Primal 


Dual 


Analytical 


Error (%) 


TMm 


2.4206 


2.4206 


2.4048 


0.6569 




3.8883 


3.8883 


3.8317 


1.4758 


TMn 


3.8901 


3.8901 


3.8317 


1.5234 


TM 21 


5.2669 


5.2699 


5.1356 


2.5563 


TM 21 


5.2694 


5.2694 


5.1356 


2.6050 


# zero modes 





311 






# nonzero modes 


136 


136 







TABLE V: Numerical results of number of modes of TE and TM. 





Primal formulation 


Dual formulation 


# zero modes (TE) 




1 


# zero modes (TM) 





N F -1 


7^ nonzero modes (TE) 


N in _ N in 


N F -1 


7^ nonzero modes (TM) 




N% - (N F - 1) 



where 5 is the co differential operator (Hilbert adjoint of d). An arbitrary contractible 2D 
domain Q can be discretized by a general grid made up of a network of polygons. We will 
briefly discuss next the connection between the discrete Hodge decomposition above and the 
Euler's formula for a network of polygons (for a more details, see reference j^l). 

1. 2D TE case 

For 2D TE case, applying to the electric field intensity E (1-form), we obtain 

E 1 = d(f) + SA 2 , (36) 

where 4>° is a 0-form and A 2 is a 2-form. In Eq. (jHEJ) d(p° represents the static field and 5 A 2 
represents the dynamic field. We can trace the following correspondence between Euler's 
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formula for a network of polygons and the Hodge decomposition 



191 



N in _ N in = Np _ 1; ( 37 ) 
I I I 

E l -d<j) = 5A 2 , 

where Ny is the number of internal vertices, N l £ the number of internal edges and Np the 
number of faces of a mesh. 

2. 2D TM case 

For 2D TM case, applying (|35|) to the electric field intensity E (0-form), we obtain 

E° = 5A\ (38) 

where A 1 is a 1-form. We can trace the following correspondence between Euler's formula 
for a network of polygons and the Hodge decomposition 

Ny — = [N™ - {N F - 1)] , (39) 

I I 
E° = 5 A 1 , 

3. Zero modes and nonzero modes 

Eq. (p?7j) or Eq. (J3T?|) can be summarized as 

L t -L 2 = R X - R 2 . (40) 

For TE case, we identify 

L\ = N%, L 2 = Ny, R, = N F ,R 2 = 1, (41) 

and for TM case, we identify 

L x = Ny 1 , L 2 = 0,R 1 = N%, R 2 = (N F - 1) . (42) 

The l.h.s. of Eq. (J4T!|) corresponds to the range space of [X E ] while the r.h.s. corresponds 
to the range space of [Xh] ■ Furthermore, the L 2 corresponds to the null space of [X E ] while 
R 2 corresponds to the null space of [Xh] ■ These results are summarized in Table VI. 
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TABLE VI: Null spaces and range spaces of [Xe] and [Xh] 





\x E ] 

L -C' J 




Dim(Null space) (TE) 




1 


Dim(Null space) (TM) 





N F -1 


Dim(Range space) (TE) 


N in _ N in 


N F -1 


Dim(Range space) (TM) 




Nf - (N F - 1) 



Table VI exactly matches Table V from numerical results. The DoFs of system matrices 
[Xe] and [Xh] equal the total number of modes of primal formulation and dual formulation, 
respectively. Furthermore, the DoFs in the null space of [Xe] and [Xh] equal the number of 
zero modes of primal formulation and dual formulation, respectively. Finally, the DoFs in 
the range space of [Xe] and [Xh] equal the number of nonzero (dynamic) modes of primal 
formulation and dual formulation, respectively. Note that in the case of 2D TE modes (the 
electric field intensity E is a 1-form interpolated by edge elements), it is a well known fact 
that the dimension of the null space (# zero modes ) of [Xe] is equal to the number of 
internal nodes SHE- 

From Eq. (pfUj) (Euler's formula for a network of polygons) it can be concluded that 
the dimension of range space of [Xe] equals the dimension of range space of [Xh], as a 
fundamental property of discrete Maxwell equations [^. 



D. Polygonal cavity 

A 2D cavity of arbitrary shape can be approximated by a polygon as the boundary P|. 
Table VII and Table VIII present the results for TE modes and TM modes of a polygonal 
cavity (Fig. 2). These results corroborate the conclusions summarized by Table V and 
Table VI. Moreover, both systems matrices [Xe] and [Xh] are finite approximation of the 
corresponding infinite system. If we use same mesh and same basis functions, that is, same 
basic matrices [d C uw], [^*«rzL K^- 1 ] an d [*<=], the dynamic physical structure encoded by 
system matrices [Xe] and [Xh] will be identical. Furthermore, if we use same linear solver, 
the solutions of both formulations will give identical nonzero modes (dynamic solutions) up 
to round off errors (see Table VII and VIII) . 
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TABLE VII: TE modes (the angular frequencies of the 5 lowest nonzero modes) of a polygonal 
cavity. 



Mode No. (TE) 


Primal formulations 


Dual formulation 


1 


2.57359064243139 


2.57359064243165 


2 


3.28134124800976 


3.28134124800987 


3 


4.32578591632893 


4.32578591632896 


4 


5.17188723866480 


5.17188723866481 


5 


5.94586993156365 


5.94586993156362 


# zero modes 


73 


1 


# nonzero modes 


175 


175 



TABLE VIII: TM modes (the angular frequencies of the 5 lowest nonzero modes) of a polygonal 
cavity. 



Mode No. (TM) 


Primal formulations 


Dual formulation 


1 


4.06172573841605 


4.06172573841600 


2 


6.20284873300873 


6.20284873300876 


3 


6.85765079948016 


6.85765079948015 


4 


8.31632816148913 


8.31632816148915 


5 


9.05550834626485 


9.05550834626483 


# zero modes 





175 


# nonzero modes 


73 


73 



VI. CONCLUDING REMARKS 

Based on Galerkin duality and discrete Hodge operators, we construct two system matri- 
ces, [X E ] (primal formulation) and [X H ] (dual formulation) that discretize the wave equa- 
tions. It can be shown that the primal formulation recovers conventional (edge-element) 
FEM and suggests a geometric foundation for it. On the other hand, the dual formulation 
suggests a new (dual) type of FEM. Although both formulations give identical physical solu- 
tions, the null spaces are different. The Hodge decomposition of the DoFs can be associated 



13 



with Euler's formula for a network of polygons in 2D or polyhedra in 3D. 
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APPENDIX A: STIFFNESS MATRIX: GEOMETRIC VIEWPOINT 

Using 3D tetrahedral and cubic elements, respectively, and assuming that the perme- 
ability n is constant within each element, we will show that stiffness matrix [S] equals the 
multiplication of incidences and Hodge matrices 



[S] = Kuril I***" 1 ] [ d curi] ■ 



(Al) 



1. 



Tetrahedral elements 



From the DoFs for the tetrahedral element (Fig. 3) 



B 



b 



1,2,3 u l,3,4 u l,4,2 "2,4,3 



(A2) 



E 



e l,2 e l,3 e l,4 e 2,3 e 4,2 e 3,4 



(A3) 



we can construct the incidence matrices [d cur i\ and [d* url ] 



1-10 10 



1-10 1 




(A4) 



-10 10 10 



-1 -1 -1 



[d*curl\ - [dcurl] ■ 



(A5) 
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In the above, the superscript t stands for transposition. Using the vector calculus proxies of 
3D Whitney 2-form, the Hodge matrix [^-i] can be calculated as 



L**» _1 J{(i,j,fc),(iJ,fc)} 

Let 



/ -^li,k ■ v?h,~k dV - ( A6 ) 



[G] = [d* curl ] [Hp-i] [d curl ] , (A7) 
which is a 6 x 6 matrix. The entry of stiffness matrix [S] can be computed as 

= - (2^0 x ^Ci) • (^(l x V^) ■ (A8) 

By comparing each term of matrix (|A7|) with the corresponding term of matrix ([ A8|) . such 
as [G] 12 

[ G ]i2 = " t^ln + K-^si + tV^ia " [V 1 ] 32 ( A9 ) 

= - (2^Ci x V^) • (2^Ct x V^) , (MO) 



and [S] 



12 



we obtain 



[S} 12 = - (2^Ci x ^(2) • (2^Ct x V^) , (AH) 



[S] = K arl ] K-i] [«U . (A12) 



2. Cubic elements 



Consider a cubic element given in Fig. 4, whose side length is L and whose center is at 
(x c ,y c ). From the DoFs for the cubic element 

1 1 



^1,4,3,2 ^5,6,7,8 ^2,3,7,6 5,8,4 ^1,2,6,5 ^3,4,8,7 



(A13) 



E 



1 t 



e l,2 e 4,3 e 5,6 e 8,7 e l,4 e 5,8 e 2,3 e 6,7 e l,5 e 2,6 e 4,8 e 3,7 



(A14) 
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we can construct the incidence matrix [d cur i] and [d* uW ] for the cubic element 



[d, 



curl\ 



-11 1-1 

1-1-1 1 

1-1-1 1 
-11 1-1 
1-1 -11 

-11 1-1 



[d*curl\ - [dcurl] 



The edge elements A^ • for a cubic element can be written as 21 1 



(A15) 



(A16) 



5,6 



^8,7 

^2,3 

^6,7 



2.(1 



1,8 



3,7 



1 

L 3 
1 

L 3 
1 

L 3 
1 

L 3 
1 

L 3 
1 

L 3 
1 

L 3 
1 

L 3 
1 

L 3 
1 

L 3 
1 

I 3 
1 

L 3 



L 



L 



Vc+ 77 +2/ 2 c + 77 ~ 2 x 



2/c + 77 - 1/ + 77 + * x > 



-2/c + 75- + 2/ 



-z c + — + z ] X, 



Z c + ^ - 2 i ( - 



L 



x y, 
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The corresponding face elements ~$i j k,i can be constructed as 



^1,4,3,2 = ~J3 {z c + ^ - zj z, 
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The Hodge matrix can be calculated as 



K-A = 



6Lfi 
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Let c = The matrix [G] can be computed as 



(A18) 



(A19) 



(A20) 



[G] = K url ] [ V ] [dcurl] 
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(A21) 

Using the formula 

'%«>.(«)} = \\ x ■ x dy - < A22 > 

the stiffness matrix [S] can be computed as 



[S] = c 
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(A23) 
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Comparison of Eq. (jA21|) and Eq. (|A23|) gives the following identity 

[S] = K url ] [d C uri] ■ (A24) 

The above proof can be straightforwardly extended to the rectangular brick element whose 
side lengths are (L x , L y , L z ). 
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Fig.l. The mesh has 178 vertices (136 internal vertices), 447 internal edges, and 312 
triangles. 

Fig.2. The coordinates of the vertices of the polygon are 

(0,0), (1,0), (1.4, 0.4), (1.3, 1.0), (0.8, 1.2), (0.3, 0.9). The mesh has 105 vertices (73 
internal vertices), 248 internal edges, and 176 triangles. 

Fig. 3. Oriented tetrahedral element. 

Fig. 4. Oriented cubic element. 
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